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ABSTRACT 

The collapse and fragmentation of filamentary primordial gas clouds are explored using one- 
dimensional and two-dimensional hydrodynamical simulations coupled with the nonequilibrium pro- 
cesses of hydrogen molecule formation. The cloud evolution is computed from the initial central density 
n c = 10 — 10 6 cm -3 . The simulations show that depending upon the initial density, there are two occa- 
sions for the fragmentation of primordial filaments. If a filament has relatively low initial density such as 
n c < 10 5 cm -3 , the radial contraction is slow due to less effective H2 cooling and appreciably decelerates 
at densities higher than a critical density, where LTE populations are achieved for the rotational levels 
of H2 molecules and the cooling timescale becomes accordingly longer than the free-fall timescale. This 
filament tends to fragment into dense clumps before the central density reaches 10 8-9 cm -3 , where H2 
cooling by three-body reactions is effective and the fragment mass is more massive than some tens M Q . 
In contrast, if a filament is initially as dense as n c J> 10 5 cm~ 3 , the more effective H2 cooling with the help 

of three-body reactions allows the filament to contract up to n ~ 10 12 cm -3 . After the density reaches 
n ~ 10 12 cm" 3 , the filament becomes optically thick to H2 lines and the radial contraction subsequently 
almost stops. At this final hydrostatic stage, the fragment mass is lowered down to w 1M Q because of 
the high density of the filament. The dependence of the fragment mass upon the initial density could be 
translated into the dependence on the local amplitude of random Gaussian density fields or the epoch 
of the collapse of a parent cloud. Hence, it is predicted that the initial mass function of Population III 
stars is likely to be bimodal with peaks of « 10 2 M Q and « lM©, where the relative heights could be 
a function of the collapse epoch. Implications for the metal enrichment by Population III stars at high 
redshifts and baryonic dark matter are briefly discussed. 

Subject headings: cosmology: theory — galaxies: formation — hydrodynamics — ISM: clouds — stars: 
formation 



1. INTRODUCTION 

The existence of a very first generation of objects, 
namely Population III, has been originally postulated by 
the presence of a noticeable amount of heavy elements 
in Population II stars (see, e.g., Carr, Bond, & Arnett 
1984 and references therein) and has recently gained in- 
creasing importance owing to the discovery of intergalac- 
tic metals in the Lya forest (Cowie et al. 1995; Songaila 
& Cowie 1996). In those respects, Population III ob- 
jects need to be such massive stars that they can pro- 
duce metals at the end of their evolution. Recent theo- 
retical analyses on the evolution of metal-free stars pre- 
dict that the fate of the massive metal-free stars can be 
classified as follows (e.g., Heger, Woosley, & Waters 2000; 
Chiosi 2000; see also Portinari, Chiosi, & Bressan 1998 
for the effects of mass loss.): (1) A star with a mass of 
m ;> 250M Q collapses completely to a black hole (BH) 
without ejecting any heavy elements (Bond, Carr & Ar- 
nett 1984). (2) A star of 100M© £ m £ 25OM is partly 
or completely disrupted by electron-positron pair instabil- 
ity. For m ^ 150M©, the core completely disrupts, and 
the whole core involving heavy elements is injected in the 
intergalactic medium. If it is an extremely energetic event 
(a hypernova), it might lead to a gamma ray burst (GRB). 



(3) A star of 35M Q <, m <, 100A/ Q probably collapses into 
a black hole. (4) A star of 10-M© <; m <J 35M Q results in 
a type II supernova. 

Population III stars are related to various issues that 
are currently the object of considerable attention. The 
luminous Population III stars could cause the reioniza- 
tion of the universe at redshifts J> 5 (Couchman & Rees 
1986; Fukugita & Kawasaki 1994; Ostriker & Gnedin 1996; 
Gnedin & Ostriker 1997; Haiman & Loeb 1997; Miralda- 
Escude & Rees 1998; Gnedin 2000). Alternatively, moder- 
ately massive BHs as the end products of massive stars 
might coagulate into a super-massive BH, evolving to 
primordial AGNs (Larson 2000). The accreting super- 
massive BHs may be more responsible for cosmic reion- 
ization (Tegmark & Silk 1994; Sasaki & Umemura 1996; 
Haiman & Loeb 1998; Valageas & Silk 1999). In addition, 
Population III stars may play an important role in the 
early evolution of galaxies (e.g., Tegmark, Silk, & Blan- 
chard 1994; Ostriker & Gnedin 1996) or the early for- 
mation of massive BHs of ~ 10 5 Mq (Umemura, Loeb, 
&: Turner 1993). They may be responsible for the ob- 
served abundance patterns of extremely metal-deficient 
stars (McWilliam et al. 1995; Ryan, Norris, & Beers 1996; 
Audouze & Silk 1995; Shigeyama & Tsujimoto 1998). Fi- 
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nally, if a significant number of MACHOs are ancient white 
dwarfs (Mcndcz & Minniti 2000), they may stem from low- 
mass Population III stars (Carr 1994; Larson 1998; Rees 
1999). In the light of all of such possible significant conse- 
quences, the initial mass function (IMF) of Population III 
stars is an issue of undoubted importance. 

Many authors have studied the collapse of primordial 
clouds to estimate the masses of Population III stars (e.g., 
Matsuda, Sato, & Takeda 1969; Yoneyama 1972; Hutchins 
1976; Silk 1977, 1983; Carlberg 1981; Palla, Salpeter, & 
Stahler 1983; Yoshii & Saio 1986; Nishi et al. 1998). These 
studies have emphasized the importance of radiative cool- 
ing by H 2 molecules because the primordial gas is deficient 
in heavy elements, which are the most efficient coolants in 
present-day star formation. In the first collapsed objects, 
baryonic gas is heated to temperatures above 10 3 K during 
the contraction. This enhances the H 2 formation rate and 
causes the H 2 abundance to rise from its initial value of 
<~ 10~ 6 to a quasi-equilibrium value of 10~ 4 ~ 10~ 3 (e.g., 
Palla et al. 1983; Lepp & Shull 1984; Haiman, Thoul & 
Loeb 1996). Thus, H 2 molecules cool the gas to a temper- 
ature of ~ 500 K. As a result, the Jeans mass descends 
to a stellar mass and Population III stars can form in this 
way through fragmentation of the first collapsed objects. 
Although many elaborate analyses have been made, the 
estimated masses of Population III stars have not been 
well converged. Several authors have suggested that Pop- 
ulation III stars were low mass, whereas others have sug- 
gested that they were massive or very massive. This dis- 
crepancy seems to come from most studies being restricted 
to highly simplified models such as homogeneous, pressure- 
less, and/or spherical collapses. 

In the bottom-up scenarios like cold dark matter (CDM) 
models, the first collapsed pregalactic objects should form 
at redshifts of z ~ 10 — 10 2 and have masses of 10 5 ~ 
10 8 M Q (Tegmark et al. 1997). Because of the asymp- 
totic scale-invariance of the CDM density fluctuations, 
pregalactic clouds with this mass range could undergo a 
run-away collapse, forming cores of <~ 10 2 M Q (Abel et 
al. 1998; Abel, Norman, & Bryan 2000), or form mini- 
pancakes, fragmenting into pieces of 10 2 M Q — 10 3 M Q 
(Bromm, Coppi, & Larson 1999). In these calculations, 
one of the common features is the formation of a filamen- 
tary structure. Filamentary clouds are gravitationally un- 
stable and likely to fragment into dense clumps. Such 
dense clumps are expected to evolve into Population III 
stars. The physics of the fragmentation of filamentary pri- 
mordial clouds has been studied analytically (Uchara et 
al. 1996) or by a one-dimensional numerical simulation 
(Nakamura & Umcmura 1999, hereafter Paper I). These 
studies suggest that the minimum mass of Population III 
stars is of the order of Mq. However, the physical pro- 
cesses are still unclear to determine whether Population 
III stars can be far above M Q or be eventually reduced 
to a few Mq. Thus, as a further step following Paper I, 
we here perform two-dimensional hydrodynamical simu- 
lations. Attention is focused on elucidating the physical 
process of the fragmentation of primordial gas filaments to 
assess the mass of Population III stars. 

In Paper I, we pursued the radial contraction of primor- 
dial filaments and showed that the filaments continue to 
contract quasi-statically, the temperatures staying nearly 
constant at T <~ 500 K. When the cloud becomes opti- 



cally thick to the H 2 lines, the radial contraction essen- 
tially stops. Applying a linear stability analysis, the frag- 
mentation was expected to take place at that stage and the 
minimum masses of Population III stars were estimated as 
a few Mq. Thus, Population III stars were anticipated to 
be low-mass deficient compared to the present-day stars. 
In this paper, we pursue the fragmentation processes of 
the filaments by means of two-dimensional axisymmetric 
simulations. The present model is an improved version of 
Paper I. The numerical model and method are described 
in §2. Numerical results are given in §3 and §4. We show 
that filaments with low initial density can fragment into 
dense clumps before the cloud becomes optically thick to 
the H 2 lines. Then, the masses of the clumps could be 
much more massive than a few Mq. However, relatively 
dense filaments result in clumps of a few Mq. Hence, in §5, 
it is predicted that the initial mass function of Population 
III stars is likely to be bimodal. In §6, we discuss some 
implications for the first metal enrichment and baryonic 
dark matter. 

2. MODEL AND NUMERICAL METHODS 
2.1. Basic Equations 

To pursue the collapse and fragmentation of filamentary 
primordial clouds, we employ a two-dimensional hydrody- 
namical scheme. We assume that the system is axisym- 
metric and that the medium consists of ideal gas. The 
adiabatic index, 7, is taken to be 5/3 for monatomic gas 
and 7/5 for diatomic gas. We deal with the following 9 
species: e, H, H+, H", H 2 , H+, He, Hc+, and Hc++. The 
mass fraction of He is taken to be 0.24 of the total mass. 

The basic equations are then described as 
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where p is the mass density, n is the number density, v is 
the velocity, P is the gas pressure, tp is the gravitational 
potential, E is the energy per unit volume, G is the grav- 
itational constant, and k is the Boltzmann constant. The 
symbol A nct denotes the net energy-loss rate per unit vol- 
ume. The values with subscript i denote those of the i-th 
species. 

The number density of the z-th species, rij, is obtained 
by solving the following time-dependent rate equations, 
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where the reaction rate coefficients, kjk and ki mn , are given 
in Table |]. The relative abundances of hydrogen species 
are given by Xi = rii/n, where n = nn + + n n- + 
2n H2 + 2n„+. 

We take into account the following thermal processes 
by H atoms and H2 molecules: (1) H cooling by recom- 
bination, collisional ionization, and collisional excitation 
(Cen 1992), (2) H2 line cooling by rotational and vibra- 
tional transitions (see below) , (3) cooling by H2 collisional 
dissociation (Shapiro & Kang 1987), and (4) heating by 
H2 formation (Susa, Uehara, & Nishi 1996). In Paper 
I, for the H2 formation heating, we only took into ac- 
count the contributions by two-body reactions (Shapiro 
& Kang 1987). In this paper, we include the contributions 
by three-body reactions, which play an important role in 
temperature evolution after the density reaches 10 9 cm~ 3 
(see §3). Other thermal processes are negligible because 
the gas temperatures did not exceed 10 4 K in the models 
calculated in this paper. 

As shown by many authors, H2 line cooling is most ef- 
ficient in primordial gas. Therefore, careful evaluation of 
the H2 line cooling rate is necessary. We thus compute 
the H2 line cooling rate as follows. First, the level pop- 
ulations at the rotational and vibrational excitations are 
determined by using a recursion formula (Hutchins 1976; 
Palla et al. 1983). For the collisional deexcitation rates, 
we consider both H-H2 and H2-H2 collisions using the an- 
alytical fits of Hollenbach & McKee (1979) and Galli & 
Palla (1998). The 21 rotational and 3 vibrational levels 
are taken into account. Next, a photon escape probability 
method is applied for each transition (Castor 1970; Gol- 
dreich & Kwan 1974). Finally, the line cooling rate is 
computed as 



In bottom-up scenarios such as CDM models, the first 
collapsed objects at z ~ 10 — 100 are expected to have 
masses of 1O 5_8 M0. In these clouds, the gas is heated 
above T > 10 3 ~ 5 K by shock. Therefore, just after 
the shock formation, the cooling timescale is likely to be 
shorter than the free-fall one. Then, the gas is cooled 
by the H2 cooling down to the temperature at which the 
cooling timescale is comparable to the free-fall one (e.g., 
Haiman et al. 1996; Yoneyama 1972). The resultant gas 
temperature is estimated to be T ~ 10 2-3 K, depending 
weakly on the gas density and the H 2 fraction. Recent 
numerical simulations have shown that such clouds tend 
to fragment into filaments that collapse toward the major 
axes (e.g., Bromm et al. 1999; Tsuribe 2000). The collaps- 
ing filaments are expected to fragment into denser clumps 
where Population III stars will form. A model of such a 
filamentary cloud is described below. 

The model of a filamentary cloud is basically the same 
as that presented in Paper I. We consider an infinitely long 
cylindrical gas cloud that is collapsing in the radial direc- 
tion. The initial temperature and relative abundances are 
assumed to be spatially constant. At the initial state, the 
relative abundances of H~ , H^ , He + , and He ++ are set to 
zero for simplicity. We do not consider the effect of dark 
matter because after virialization of the parent system, 
the local density of baryonic gas is likely to become higher 
than the background dark matter density owing to radia- 
tive cooling (e.g., Cen & Ostriker 1992a, 1992b; Umemura 
1993). The density is assumed to be uniform along the 
cylinder axis and its radial distribution is expressed as 
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where n^i is the H2 level population at level i, Aij is the 
Einstein A-coefficient, and /iz/y is the energy difference be- 
tween levels i and j. The photon escape probability /3y is 
defined as 

N - iz^zZk) , ( io) 
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where T,_>j is the optical depth at the transition i — > j, 
Tout is the radius of the cloud surface, is the absorp- 

tion coefficient of the transition i — > j, Bij and Bji are the 
Einstein B-coefficients, and Az/y(oc T 1 / 2 ) is the thermal 
Doppler width of the transition line i — » j. Our cooling 
function almost coincides with that of Galli & Palla (1998) 
as long as the cloud is optically thin to the H2 lines. How- 
ever, once the cloud becomes opaque to the cooling radi- 
ation, our cooling rate significantly deviates from that of 
Galli & Palla (1998). 

2.2. Model of Filamentary Primordial Clouds 



where Rq = ^/2/fcTo/ (^Gpop) is the effective radius, po is 
the central mass density, T is the initial gas temperature, 
p is the mean molecular weight, and / is the ratio of the 
gravitational force to the pressure force. When / = 1, the 
density distribution coincides with that of an isothermal 
filament in equilibrium (Stodolkiewicz 1963). In this pa- 
per, we restrict the parameter / to / > 1 since we are 
interested in the evolution of collapsing clouds. 
The radial infall velocity is given by 



v r 



(14) 



which has a qualitatively similar distribution to that of 
a self-similar solution of a collapsing filament (see Figure 
6 of Paper I). Here, i>o is constant and is set to an initial 
sound speed (c s ) for all the models calculated in this paper. 
We also calculated the evolution of the models with dif- 
ferent velocity distributions, e.g., v r = — vq sin(7rr/i? max ) 
and confirmed that the numerical results are not sensitive 
to the assumed velocity profiles as long as vq <; a few x c s . 
This is because even if the initial cloud is static, the ra- 
dial contraction is immediately accelerated owing to the 
H 2 cooling. 

For the above model, the mass per unit length (line 
mass) is given by 
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= l^x 10»M aP c-. (0 (^) . (15) 

Note that the line mass of the equilibrium (/ = 1) fila- 
ment depends only on the temperature. When the fila- 
mentary cloud forms through the gravitational fragmenta- 
tion of a parent sheetlike cloud, its line mass is expected 
to be nearly twice that of the equilibrium cloud (see §2.2 
of Paper I). We thus adopt / ~ 1.5 for typical models. 

This model is then specified by the five parameters: 
tiq = po/Mj la, f, x e , and ih 2 - The abundances of x e , 
in, and i H + are determined by the conservation of mass 
and charge. The initial density of the filament depends 
sensitively on the properties of the parent cloud such as 
the mass, collapse redshift, and spin parameter (see the 
Appendix). Thus, we take a wide range of the initial den- 
sity, 10 cm" 3 < no < 10 6 can . Previous numerical sim- 
ulations (e.g., Haiman, Rees, & Loeb 1996; Abel et al. 
1998; Bromm et al. 1999) have shown that the temper- 
atures and H2 relative abundances reach T ~ 300 — 500 
K and Xn 2 ~ 10~ 3 — 10~ 4 , respectively, when the first 
collapsed objects with masses of 10 5_8 M© are virialized. 
Thus, the initial temperatures and H2 abundances are set 
to T = 300 - 500 K and x Ha ~ 10~ 3 - 10~ 4 , respectively. 
The electron abundance is set to 5 x 10 -5 for all models. 
As mentioned in Paper I, the numerical results do not de- 
pend sensitively upon x e as far as x e J> 10 , because free 
electrons can quickly recombine to a level of x e < 10~ 7 in 
the course of the collapse. 

2.3. Density Fluctuations and Initial Parameters 

The initial model of the computation is constructed by 
imposing linear density fluctuations on the model filament 
described in the above subsection. The power spectrum 
of the density fluctuations is assumed to be a power law 
distribution of P(k) = Ak v on a wavenumber space, where 
k is the wavenumber of the density fluctuation in the in- 
direction and v is the power index. The phase of the fluc- 
tuations has a random distribution in the range of to 2ir. 
In this paper, the power index v is taken to be v — —2, — 1, 
and 0. (Note that v = — 1 provides scale-invariant fluctua- 
tions in the present geometry.) The amplitude is specified 
by 5 = [P^k] 1 / 2 at k = 2n/R . The simulations were 
done in the range of 0.1 > 5 > 0.001. 

Table summarizes the initial parameters of the numer- 
ical models calculated in this paper. 

2.4. Numerical Methods 

The hydrodynamic equations (Q) through (^) are solved 
numerically using a second-order upwind scheme based on 
the method by Nobuta & Hanawa (1999). See Nobuta & 
Hanawa (1999) for more details and the test of the code. 
The rate equations (||) are solved numerically with the 
LSODAR (Livemore Solver for Ordinary Differential equa- 
tions with Automatic method switching for stiff and non- 
stiff problems) coded by L. Petzold and A. Hindmarsh. 

The computations are performed on a cylindrical do- 
main (0 < r < -Rmax and < z < Z max ) with a fixed 
boundary condition at r = i? max and periodic boundary 
conditions at z — and Z max , where i? max and Z max are 
the maximum r- and z-coordinates on the computational 
domain, respectively. Uniform grids are employed for all 



runs. In all the models we take i? max to be more than 
three times the effective radius Ro- The effect of the fixed 
boundary is very small. This is because the density is 
much lower near the outer boundary (r ~ i? m ax) than at 
the center, i.e., p <; 10 _3 p c for all runs. In fact, we have 
calculated the cases with larger i? max and have confirmed 
that the numerical results are not affected. 

3. RADIAL CONTRACTION OF FILAMENTARY CLOUDS 

Before examining the fragmentation processes of fila- 
mentary clouds, we pursue the evolution of the filamen- 
tary clouds with no density fluctuations, i.e., their radial 
contraction. These calculations are the extension of Paper 
I. (In Paper I, we explored the evolution of the filaments in 
the initial density range of 10 2 cm~ 3 < uq < 10 4 cm~ 3 . In 
this paper, we extend the initial density range to 10 cm 
< no < 10 6 cm -3 .) In the following we focus on the evo- 
lution of physical quantities on the z-axis. See Paper I for 
more details, e.g., the radial distributions of the physical 
quantities. As shown below, the evolution of the filaments 
can be classified into two types, depending mainly on the 
initial density. 

3.1. A Low-Density Filament: A Quasi- Statically 
Contracting Filament (Model Ala) 

As a typical example of low-density filaments, we show 
the evolution of model Ala which has the initial parame- 
ters of n = 10 cm" 3 , T = 400 K, and / = 1.5. The elec- 
tron and H2 number fractions are initially set to 5 x 10~ 5 
and 10~ 4 , respectively. Figure la shows the evolution of 
the temperature and the H2 number fraction as a function 
of the central density. The central density monotonously 
increases with time. Therefore, the abscissa corresponds to 
the evolution time. For comparison, we pursued the evo- 
lution of model Ala in which the three-body reactions are 
not taken into account. The evolutional paths are shown 
by dotted lines in Figure 1. 

Figure lb shows the evolution of three characteristic 
timescales as a function of the central density. The solid, 
dashed, and dashed-dotted lines denote the contraction 
time, the cooling time, and the fragmentation time, re- 
spectively. The dotted lines denote the evolution of the 
temperature and contraction time for the model without 
the three-body reactions. The contraction and cooling 
times are defined as td yn = p/f> and i coo i = 3nfcT/(2A nct ), 
respectively. The fragmentation time is defined to be the 
inverse of the growth rate of the fastest-growing linear per- 
turbation, iftag = 2.07/ y / 2nGp c (eq. [38] of Nakamura, 
Hanawa, & Nakano 1993). It should be noted that a fila- 
ment does not undergo fragmentation in one-dimensional 
calculations, therefore if rag should be regarded as a mea- 
sure of the free-fall time. 

During the contraction, the temperature stays nearly 
constant at T ~ 300 — 500 K because of the H2 line cool- 
ing. The H2 number fraction also stays nearly constant 
at Xb_ 2 — 10~ 4 ~ 10~ 3 until the density reaches n c ~ 10 8 
cm~ 3 . After that, the three-body reactions for the H2 
formation become dominant. Accordingly, the H2 num- 
ber fraction steeply rises around n c ~ 10 s cm -3 and al- 
most all the hydrogen atoms become into H2 molecules by 
the stage at which the density reaches n c ~ 10 11 cm~ 3 . 
In contrast, for the model without the three-body reac- 
tions, the H2 number fraction stays nearly constant at 
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xa 2 = 10 -4 ~ 10 -3 and the contraction stops at the stage 
at which the density reaches n c ~ 10 9 cm -3 . 

The heating rate by H2 formation also becomes appre- 
ciable when the three-body reactions become dominant. 
Therefore, the temperature rises slowly after the density 
reaches n c ~ 10 9 cm -3 . Note that the H2 line cooling rate 
is always the most efficient during the contraction. 

During the contraction, the contraction time nearly co- 
incides with the cooling time until the central density 
reaches 2 x 10 11 cm -3 . This indicates that the cloud col- 
lapses in a cooling time. After the central density reaches 
fO 4 cm -3 , the contraction proceeds quasi-statically be- 
cause the contraction time becomes longer than the frag- 
mentation time. In other words, if there are fluctuations, 
the cloud is expected to become unstable to fragmentation 
after the density reaches 10 4 cm -3 . Such evolution is ex- 
plained as follows. A density of 10 4 cm -3 is almost compa- 
rable to a critical density of H2, n cr ~ 10 3-4 cm -3 , beyond 
which LTE populations are achieved for the rotational lev- 
els of hydrogen molecules. If the density is less than n cr , 
then the H2 line cooling rate is nearly proportional to the 
square of the density, Aotn 2 , and the cooling time is thus 
inversely proportional to the density, t coo \ cx n . On the 
other hand, if the density is greater than n cr , then the 
H2 line cooling rate is nearly proportional to the density, 
A cx n and accordingly the cooling time is independent of 
the density. Therefore, after the density reaches n cr , the 
cooling time becomes longer than the fragmentation time 
(tfrag oc n -1 / 2 ) when the temperature is nearly constant. 
(In fact, the cooling time is almost constant during 10 
cm -3 <; no < 10 7 cm -3 .) After the central density reaches 

10 s cm -3 , the H2 abundance steeply increases due to the 
three-body reactions and the radial contraction thus ac- 
celerates again owing to the enhanced H2 line cooling. (In 
contrast, for the model without three-body reactions, the 
radial contraction stops when the central density reaches 
~ 10 9 cm -3 because of less effective H2 cooling.) When the 
central density reaches 10 11-12 cm -3 , the cloud becomes 
optically thick to almost all the H2 lines which signifi- 
cantly contribute to the total cooling rate. Consequently, 
the contraction essentially stops at that stage. 

Such evolution is also related to the dynamical stabil- 
ity of self-gravitating clouds. Assuming the polytropic 
relation of the equation of state (P oc p 7 ), a hydro- 
static cylindrical cloud is stable to radial contraction when 
7 > 7cr = 1. In the present model, the temperature is 
nearly constant during the contraction. In other words, 
the cloud is in a marginally stable state. Thus, once the 
H2 line cooling becomes less effective, the radial contrac- 
tion easily stops. This is essentially different from that of 
a spherical cloud. A spherical polytropic gas cloud can 
continue to collapse as long as 7 < 7 cr = 4/3 even if the 
cooling rate becomes less effective (see §4.1 and Omukai 
& Nishi 1998). 

3.2. A Dense Filament: A Dynamically Collapsing 
Filament (Model C6a) 

As a typical example of dense filaments, we show the 
evolution of model C6a which has the initial parameters 
of n = 10 6 cm" 3 , T = 400 K, and / = 4. The elec- 
tron and H2 number fractions are initially set to 5 x 10 -5 
and 10 -4 , respectively. Figures 2a and 2b are the same 



as those of Figures la and lb, respectively, but for model 
C6a. 

The evolution is qualitatively similar to that of model 
Ala, although the temperatures are about two times 
higher during the contraction. Since the H2 formation 
rates by the three-body reactions are inversely propor- 
tional to the temperature, the three-body reactions be- 
come more effective at the later stages than for model 
Ala (n > 10 9 cm -3 ). The contraction time does not be- 
come longer than the fragmentation time until the cloud 
becomes optically thick to the H2 lines (n c ~ 4 x 10 12 
cm -3 ). In other words, the contraction proceeds dynami- 
cally. This is due to the effective cooling by the three-body 
reactions. Therefore, the fragmentation is not expected to 
take place until the cloud becomes optically thick to the 
H 2 lines (n c ~ 10 12-13 cm -3 ). Note that, in the model 
without the three-body reactions, the contraction time be- 
comes longer than the fragmentation time after the central 
density reaches ~ 10 10 cm -3 , and the contraction stops 
when the central density reaches ~ 10 11 cm -3 , i.e., before 
the cloud becomes optically thick to the H2 lines. 

3.3. Summary of One-Dimensional Simulations 

We pursued the evolution for the other model parame- 
ters tabulated in Table g. We found that the evolution of 
all the models falls into either of the two types described 
above, depending mainly on the initial density. When the 
initial density is lower than no <; 10 5 cm -3 , the contrac- 
tion time becomes longer than the fragmentation time be- 
fore the three-body reactions become effective (n ~ 10 8-9 
cm -3 ). However, the contraction time does not exceed 
the fragmentation time until the density becomes greater 
than n cr , beyond which LTE populations are achieved for 
the rotational levels of H2 molecules. On the other hand, 
when the initial density is greater than no > 10 5 cm -3 
and / ;> 3, the contraction proceeds dynamically until the 
cloud becomes optically thick to the H2 lines (n c ~ 10 12 
cm -3 ). For all the models, the radial contraction essen- 
tially stops at the stage at which the cloud becomes op- 
tically thick to the H2 lines which significantly contribute 
to the total cooling rate. 

From the numerical results of the one-dimensional sim- 
ulations, the filament is likely to fragment into pieces dur- 
ing the stages at which the central density is greater than 
~ 10 4 cm -3 (p> n cr ) and is less than ~ 10 12 cm -3 because 
the radial contraction is appreciably decelerated. When 
the filament has an initial density lower than 10 5 cm -3 , 
the fragmentation is expected to take place by the stages 
at which the three-body reactions become effective. On 
the other hand, when the filament has an initial density 
greater than 10 5 cm -3 , the fragmentation is expected to 
take place at the stages at which the cloud becomes opti- 
cally thick to the H2 lines (n c ~ 10 12 cm -3 ). 

In the next section, we pursue the fragmentation pro- 
cesses with two-dimensional simulations and evaluate the 
masses of the fragments. 

4. FRAGMENTATION OF FILAMENTARY CLOUDS 

In this section, to estimate the masses of the fragments, 
we explore the evolution of collapsing filamentary clouds 
with density fluctuations by means of two-dimensional 
simulations. 
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As shown below, for several models (particularly for 
models with a large /), radial contraction proceeds to a 
great degree before the density fluctuations grow nonlin- 
early. The spatial resolution thus becomes poor before 
the cloud fragments into smaller clumps. To resolve the 
fragmentation, we refine grids according to the following 
procedure: (1) First, in the linear stage of density fluctu- 
ations, the cloud evolution for a model is pursued parallel 
by a one-dimensional calculation without density fluctu- 
ations and by a two-dimensional calculation with density 
fluctuations. The grid spacings used in a one-dimensional 
calculation are taken to be ten times finer than those of a 
two-dimensional calculation. The one-dimensional calcu- 
lation monitors the radial density profiles fairly accurately, 
while the two-dimensional calculation traces the growth 
of density fluctuations in the linear stage. (2) When the 
mean density in the z-axis reaches a reference density p Te { , 
which ensures the linear stage (typically p rc i = 10 3-4 /9o), 
the density fluctuations in the two-dimensional calculation 
are Fourier-transformed to give the power spectrum of the 
fluctuations. (3) Then, the grids in the two-dimensional 
calculation are refined in the region of < r < 0.1i? max 
and < z < 0.1Z max and the radial density profiles by 
the one-dimensional calculation are mapped upon the re- 
fined grids. (4) Finally, the Fourier-transformed density 
fluctuations obtained above are added on the refined two- 
dimensional grids. Then, the evolution of that model is 
pursued on the new computational domain. The num- 
ber of grids is 1024 2 for usual cases and 2048 2 for high- 
resolution cases. 

4.1. A Low-Density Filament (Model AJ^a) 

In this subsection, we show the evolution of model A4a 
as a typical example of less dense filaments. Figure 3 shows 
the cross-sections of the cloud at four different stages. 
This model has the initial parameters of no = 10 4 cm -3 , 
To = 400 K, and / = 1.5. At the initial state, scale- 
invariant density fluctuations [y = — 1) with an amplitude 
of Sp/p = 0.1 were added. The grid number was set to 
2048 2 . 

At the early stages, the density fluctuations do not grow 
appreciably in time because the contraction time is shorter 
than the fragmentation time [panel (b) of Fig. 3] . When 
the mean density in the z-axis exceeds ~ 7x 10 6 cm' 3 , the 
contraction time becomes longer than the fragmentation 
time and the density fluctuations begin to grow nonlin- 
early. As a result, the filamentary cloud fragments into 
denser clumps by the stages at which the mean density 
in the z-axis reaches 4 x 10 7 cm -3 [panels (c) and (d) of 
Fig. 3]. The masses of the clumps are estimated to be 
~ 250A/©. The mean separation of the clumps is nearly 
equal to 0.24 pc which is comparable to the wavelength of 
the fastest-growing perturbation at the stage at which the 
mean density in the z-axis reaches 10 7 cm' 3 . 

To see the structures of the clumps more quantitatively, 
we show the density and temperature distributions in the 
z-axis in Figure 4. At the early stages, the temperature 
stays nearly constant at T ~ 500K. After the mean den- 
sity in the z-axis reaches 10 7 cm -3 , the density fluctu- 
ations grow nonlinearly, and dense prolate clumps form. 
As the collapse proceeds, the central region of the clump 
becomes spherical. In the clumps, the contraction is ac- 
celerated, and the temperature rises slowly. This accelera- 



tion is related to the dynamical stability of self-gravitating 
clouds. A cylindrical poly tropic (P oc p 1 ) cloud is stable 
to radial contraction when 7 > 7 cr = 1. Therefore, the 
primordial filament collapses quasi-statically because the 
effective 7 is slightly greater than 7 cr (see Fig. 1). On the 
other hand, for a spherical cloud, the critical value of 7 is 
equal to 7 cr = 4/3. Accordingly, once the fragmentation 
takes place, the clumps become unstable to dynamical con- 
traction, resulting in a temperature rise. Such evolution 
is similar to that of the spherical collapse of primordial 
clouds (e.g., Omukai et al. 1998; Omukai & Nishi 1998). 

The evolutions of other models with no <; 10 5 ' 6 cm -3 
or / <J 2 are qualitatively similar to that of this model 
(e.g., models Ala - A6a, Bla - B6a, Cla - C3a). The 
fragmentation can take place before the stages at which 
the three-body reactions become efficient and the radial 
contraction is reaccelerated. 

4.2. A Dense Filament (Model C6a) 

In this subsection, we show the evolution of model C6a 
as a typical example of dense filaments. Figure 5 shows 
the density and temperature distributions in the z-axis at 
three different stages. This model has the initial param- 
eters of n = 10 6 cm" 3 , T = 400 K, and / = 4. At the 
initial state, scale- invariant density fluctuations (y = — 1) 
with an amplitude of Sp/p = 0.1 were added. As expected 
from the numerical results of the one-dimensional simu- 
lations, the radial contraction proceeds dynamically un- 
til the mean density in the z-axis reaches 4 x 10 12 cm' 3 . 
When the mean density in the z-axis exceeds 4 x 10 12 
cm -3 , the cloud becomes optically thick to the H2 lines, 
and the density fluctuations then begin to grow nonlin- 
early. In this way, the cloud fragments into clumps. The 
fragment mass is reduced down to 1 — 2Mq owing to the 
high density of the filament. 

The evolutions of other models with no p> 10 5 ' 6 cm' 3 
and / ;> 2 are similar to that of this model (e.g., models 
C4a — C6a). In those models, the fragment masses take 
their minimum values of 1 — 2M© because the radial con- 
traction proceeds dynamically until the stages at which 
the H2 lines become opaque. 

4.3. Typical Masses of the Fragments 

Figure 6 shows the distribution of the averaged frag- 
ment mass derived from the two-dimensional simulations 
for the models with v = — 1 and Sp/p — 0.1. The ab- 
scissa and ordinate denote the initial central density and 
the parameter /, respectively. Since for all the models cal- 
culated in Figure 6, the initial temperatures are taken to 
be constant at To = 400 K, the ordinate specifies the initial 
line mass. The solid lines denote the contours of the av- 
eraged fragment mass. As discussed in the Appendix, the 
primordial filaments are expected to form by cosmological 
pancake collapse and fragmentation. For comparison, the 
line masses of such filaments with I = l cq , I = 1.5Z eq , and 
I = 2l eq are shown by dashed lines in Figure 6 (see eq.]l5]|). 
Here, l cq is the line mass of the filament in hydrostatic 
equilibrium and is determined by the gas temperatures at 
which the cooling time balances with the fragmentation 
time of the pancaking disk, where the fragmentation time 
is defined as the inverse of the growth rate of the fastest- 
growing linear perturbation (Larson 1985). 
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The averaged fragment mass depends on the initial 
parameters. For larger / and/or higher initial density, 
the fragment mass is lower. As expected from the one- 
dimensional simulations, the fragmentation takes place 
during the stages at which 10 4 cm -3 <, n c <J 10 12 cm -3 
because the radial contraction proceeds quasi-statically. 
Then, the maximum and minimum masses are estimated 
as 10 3 M Q and 1^2 M Q , respectively. It is worth noting 
that these two masses are related to the microphysics of H2 
molecules. The former corresponds to the Jeans mass at 
the stage at which the density reaches the critical density 
of H 2 molecules, while the latter corresponds to the Jeans 
mass at the stage at which the cloud becomes opaque to 
the H2 lines. 

There is a steep boundary in the fragment mass at 
n ~ 10 5-6 cm -3 and / ;> 3. For the models with n ;> 10 5 
cm -3 , the fragment masses take their minimums at 1 ~ 2 
Mq. On the other hand, for the models with no <, 10 5 
cm -3 , they are greater than ~ 10 2 M©. This sensitiv- 
ity in the fragment mass comes from the rapid increase in 
H2 abundance due to the three-body reactions. As shown 
in Figure 1, when the three-body reactions become effec- 
tive (n 10 8 cm -3 ), the radial contraction accelerates 
again because of the enhanced H2 line cooling. For mod- 
els with a smaller initial density (no <J 10 5 cm -3 ), linear 
density fluctuations can grow nonlinearly before the three- 
body reactions become dominant at n <; 10 8-9 cm -3 . On 
the other hand, for models with a denser initial density 
(no ;> 10 5-6 cm -3 ), the contraction time does not exceed 
the fragmentation time until the cloud becomes optically 
thick to the H2 lines, n ~ 10 12 cm" 3 . 

Although the fragment mass also depends on the initial 
temperature, the effect of the initial temperature is iden- 
tical to that of parameter /. This is because changing To 
and/or / corresponds to the change of the line mass (see 
equation |fl5|| ). During the contraction, the cloud temper- 
ature is determined by the balance between the heating 
and cooling rates. In our model, the main heating source 
is the compressional heating by gravitational contraction, 
while the main cooling source is the H2 line transitions. 
Even if we take higher or lower initial temperature, the 
cloud temperature settles immediately to an equilibrium 
value at which the heating rate is equal to the cooling rate. 
Furthermore, the equilibrium temperature depends only 
weakly on density (T oq - 300 - 800 K for n = 10 - 10 6 
cm -3 ). Therefore, the ordinate of Figure 6 can be replaced 
by To if a constant / is adopted. 

The fragment mass depends on the amplitude and power 
index of the density fluctuations. In Figure 7, the results 
with S = 0.01 are shown. This figure shows that the de- 
pendence on the fluctuation amplitude is quite weak and 
that the fragment mass is not basically altered. Chang- 
ing power index v leads to a change of amplitude for the 
most unstable mode. It is found that the fragment mass 
is insensitive to the choice of v in a range of —2 < v < 0. 

When the initial H2 abundance is as high as 10 -3 , the 
fragment mass is reduced by a few tens % because of the 
lower temperatures although the minimum values of the 
fragment mass do not change. 



5. IMPLICATIONS FOR THE IMF OF POPULATION III 
STARS 



As shown in the previous section, the primordial fila- 
ments fragment into dense clumps whose masses are in the 
range of 1M© M <; 10 3 M©. The masses of the clumps 
depend on the initial model parameters, particularly the 
initial density. The initial densities of the filaments are re- 
lated to the initial conditions of the parent clouds. In the 
Appendix, based on a CDM cosmology, we considered the 
formation processes of the filament s an d estimated plau- 
sible initial conditions. From eq. (A4) of the Appendix, 
the initial densities of the filaments are estimated as 10 3 
cm -3 <J no 10 7 cm -3 for ler density fluctuations with 
masses of M ~ 10 6-8 Af© (10 4 cm -3 < n < 10 8 cm -3 for 
3cr density fluctuations). Then, the clump masses are eval- 
uated as 1M Q <; M <; 500M© for la density fluctuations 
(1M© <, M <; 25OM for 3a density fluctuations). Here, 
we assumed that the radius of the parent disk ranges from 
0.01r v i r to 0.1r v ; r , where r vlr denotes the virial radius of 



the parent cloud (see eq. |A1|). The minimum radius cor- 
responds to that of a rotationally supported disk with a 
spin parameter of A = 0.05. The maximum radius is taken 
from the numerical results by Bromm et al. (1999) who 
followed the collapse of 3cr 'top- hat' density fluctuations 
with a mass of ~ 10 6 Mq. Their numerical simulations in- 
dicate that by the epoch of filament formation, the radius 
of the disk shrinks to ~ 0.1r v i r , i.e., before the disk con- 
tracts to form a rotationally supported disk, fragmentation 
takes place. 

The dense clumps are expected to be the sites of Popu- 
lation III star formation. Recently, Larson (2000) argues 
that these clumps are not likely to fragment into many 
lower-mass objects because their masses are nearly com- 
parable to the Jeans mass at the epoch of fragmentation. 
Actually, numerical simulations (e.g., Bromm et al. 1999) 
have shown that the fragmentation of collapsing Jeans- 
mass clumps is likely to be limited to the formation of 
binary or small multiple systems. 

In the primordial gas, most of the parent clump mass is 
expected to accrete onto the subclumps that will evolve 
into Population III stars because metal-free (dust-free) 
gas is impervious to strong radiation pressure (Omukai 
1999, private communication). Thus, the masses of Pop- 
ulation III stars are anticipated to become comparable to 
the masses of the clumps. 

As mentioned in the previous section, the dependence 
of the clump mass on the initial density exhibits a step 
around n ~ 10 5 cm -3 . Then, the IMF of Population III 
stars is likely to be low-mass deficient and double-peaked 



at m p i = 1 — 2Mq and m 



p2 



a few x 10 - 10 2 Af, 



©• 



[The first peak is consistent with the estimates by Ue- 
hara et al. (1996) and Nakamura & Umemura (1999). 
The clumps around the second peak have similar masses 
to those obtained by Abel et al. (1998) and Bromm et al. 
(1999).] The masses of the clumps probably increase by 
merging with themselves. The resultant mass spectrum 
could have two power-law-like components with different 
peaks of m p i ~ 1 — 2M© and m P 2 ~ a few x 10 — 10 2 M©. 
It should be noted that the relative height of the first peak 
probably descends with time compared to that of the sec- 
ond peak, because the initial densities of the filaments de- 
crease with time. In other words, at higher redshifts or 
for higher a density fluctuations, the contribution of the 
lower- mass component is more significant in the IMF. 
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6. METAL ENRICHMENT 

As discussed in §|], Population III stars are expected 
to be low-mass deficient compared to present-day stars 
and their IMF is likely to be bimodal with peaks of 
m p i ~ 1 — 2Mq and m P 2 ~ a few x 10— 10 2 M Q . In Figure 
8, we show a schematic IMF of Population III stars. 

If each component of the IMF is approximated by a 
simple power law with a sharp cut-off, then the IMF of 
Population III stars is expressed as 



dN _ d{N low + N high ) 



d\ogm 



where 



dN, 



low 



dlogm 
dN high 
dlogm 



dlogm 



Am a for m > m p i 
for m < m p i 

Bm~P for m > m p 2 
for m < m p 2 



(16) 

(17) 
(18) 



and 



where the value of m P 2 is somewhat arbitrary because it 
depends on the initial densities of the filaments. The rel- 
ative heights of two peaks are also related to the power 
spectrum of the cosmological density fluctuations. If the 
power indexes of the IMF are greater than unity (a and 
[3 > 1), then the numerical constants A and B are approx- 
imated as 

A ~ (a - l^eAftota,™^" 1 , (19) 

B~ (/3 -1)(1- K)eM tota im^2 1 , (20) 

respectively, where M to tai is the parent cloud mass, e is the 
star-formation efficiency in the cloud, and k is the ratio of 
the mass contained in the lower-mass component to the 
total stellar mass. 

According to recent theoretical studies on the evo- 
lution of metal-free stars, the massive metal-free stars 
with masses of (1) 10Af o <; M star < 35M@ and (2) 
150Mq <; M s tar <, 250Mq can enrich the intergalactic 
medium through supernova explosions (e.g., Heger et al. 
2000). Roughly speaking, the former star ejects about 
10% of its total mass (~ 1M®) as heavy elements through 
a supernova, while the latter ejects about 50% (~ 75M ). 

Then, the metallicity produced by Pop ulation III stars 
can be estimated from the IMF (eq. p6[ ) as functions of 
e, k, a, and (3. When the power indexes are tentatively 
assumed to be equal to each other (a = (3), the metallicity 
is evaluated as 



3.9 x ViT z Z ( . 



(£) 



2.2 x lO -2 ^ 



1 - K 
0.5 



10" 



for a — (3 — 1.35, 



(21) 



Z . 



3.3 x 1O-Z (JL) 



9.3 x lO" 3 ^ 



1 



Go- 2 ) 



0.5 

for a — f3 = 3 



(22) 



where we took the second peak of m P 2 = 50Mq. When 
the contribution from the high-mass component is more 



significant, the metallicity is larger. If the star formation 
efficiency is as high as the present-day value of e ~ 10 -2 
and the power indexes are around a — (3 ~ 1.35 — 3, then 
the metallicity is estimated as 10~ 3 — 10 _2 Z©, which is 
consistent with the metallicity Z <J 10~ 2 Z & observed in 
Lya forest clouds with A/hi > 3 x 10 14 cm~ 3 by Cowie et 
al. (1995) and Cowie & Songaila (1998). Thus, the heavy 
elements by the first enrichment may be responsible for the 
metallicity in the intergalactic medium at high redshifts. 
This might be testable with observations of the metallicity 
and abundance ratios of heavy elements in those objects. 

Furthermore, the first enrichment by Population III 
stars might play a significant role in the early evolution 
of galaxies or abundances of the intergalactic medium ob- 
served by X-ray (e.g., Zepf & Silk 1996; Larson 1998). 



7. CONCLUSIONS 

We have explored the collapse and fragmentation of 
filamentary primordial gas clouds numerically, including 
the nonequilibrium processes for hydrogen molecule forma- 
tion. The simulations have shown that, depending upon 
the initial density, the evolution of the filaments can be 
classified into two types. If a filament has relatively low 
initial density such as n c < 10 5 cm -3 , the radial con- 
traction is slow due to less effective H2 cooling, and it 
appreciably decelerates at density higher than a critical 
density, where LTE populations are achieved for the ro- 
tational levels of H2 molecules and the cooling timescale 
becomes accordingly longer than the free-fall timescale. 
This filament tends to fragment into dense clumps before 
the central density reaches 10 8 " 9 cm~ 3 where H2 cooling 
by three-body reaction is effective, and the clump mass is 
more massive than some tens Mq. In contrast, if a fila- 
ment is initially as dense as n c ;> 10 5 cm -3 , more effective 
H2 cooling with the help of three-body reaction allows the 
filament to contract up to n ~ 10 12 cm -3 , for which the 
filament becomes optically thick to H2 lines and then the 
radial contraction almost stops. At this final hydrostatic 
stage, the clump mass is lowered down to s=s 1 — 2M Q be- 
cause of the high density of the filament. The dependence 
of clump mass upon the initial density could be trans- 
lated into the dependence of the local amplitude of ran- 
dom Gaussian density fields or the epoch of collapse of a 
parent cloud. Hence, the distribution of the clump mass 
predicts that the initial mass function of Population III 
stars is likely to be bimodal with peaks of « 1 — 2M Q and 
w 10 2 M©, where the relative heights could be a function 
of the collapse epoch. At higher redshifts or for higher 
a density fluctuations, the contribution of the lower-mass 
component is likely to be more significant in the IMF and 
the relative height of the first peak probably decreases with 
time because the initial densities of the filaments are likely 
to descend with time. 

In our model, we do not take the effects of external ra- 
diation into account. As suggested by Abel et al. (2000), 
a single star may form at the central high-density region 
of the first collapsed low-mass objects. The radiative feed- 
back from the very first stars might be significant for sub- 
sequent star formation in the surrounding medium (e.g., 
Ferrara 1998; Ciardi, Ferrara, & Abel 2000). 

We modeled the Population III IMF as a superposition 
of two power-law distributions with different peaks and es- 
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timated the metallicity produced by Population III stars 
in the first collapsed objects at high redshifts. If the star 
formation efficiency is of the same order as the present- 
day value e ~ 10~ 2 , then the metallicity is estimated as 
10~ 3 — 10~ 2 Zq. This metallicity is consistent with that 
observed in the intergalactic medium at high redshifts. 

If a significant amount of stars with masses between a 
few Mq and 8M were formed, they might have evolved 
into white dwarfs until the present epoch. The old white 
dwarfs might presently reside in the galactic halo and may 
be related to dark component observed in microlensing ex- 
periments, i.e., massive compact halo objects (MACHOs). 
If the star formation efficiency is of the same order as the 
present-day value, then the ancient white dwarfs are likely 
to contribute a few tens % of the dark mass in the galactic 
halo, which seems to be consistent with some constraints 



discussed by several authors (e.g., Chariot & Silk 1995; Al- 
cock et al. 2000; Lasserre et al. 2000; Mendez & Minniti 
2000; Hodgkin et al. 2000; Ibata et al. 2000). 
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APPENDIX 

FORMATION OF FILAMENTARY PRIMORDIAL CLOUDS 

In this appendix, we consider the formation process of filamentary primordial clouds to estimate the typical values as 
the initial conditions. 

We premise a gravitational instability scenario for the formation of cosmic structure. A cosmological density pertur- 
bation larger than the Jeans scale at the recombination epoch forms a flat pancake-like disk. This process has been 
extensively studied by many authors (e.g., Zel'dovich 1970; Sunyaev & Zel'dovich 1972; Cen & Ostrikcr 1992a, 1992b; 
Umemura 1993). Although the pancake formation was originally studied by Zel'dovich (1970) in the context of the adia- 
batic fluctuations in baryon or hot-dark-matter-dominated universes, recent numerical simulations have shown that such 
pancake structures also emerge in CDM cosmology (e.g., Cen et al. 1994). Thus, the pancakes are thought to be a 
ubiquitous feature in gravitational instability scenarios. In the following, we reconsider the pancaking of a cosmological 
density perturbation and fragmentation of the pancake into filamentary clouds. 

We first suppose a spherical top-hat overdense region in the Einstein-de Sitter Universe. The overdense region collapses 
due to self-gravity until a shock develops. After thermalization by the shock, the region forms a virialized system. The 
virial radius and temperature of the overdense region are then estimated as (e.g., Padmanabhan 1993) 

and 

where ho. 5 is the Hubble constant in units of 50 km s _1 Mpc -1 , z v ; r is the redshift epoch of virialization, Mb is the 
baryonic mass, and Oh is the baryonic density parameter. 

In a realistic situation, however, the overdense region is more or less aspherical. The deviation from spherical symmetry 
grows with time and a pancake disk consequently forms. After the shock is thermalized, the pressure force nearly balances 
the gravitational force in the vertical direction in the disk. If the radiative cooling is effective, the temperature can descend 
to a lower value than the virial temperature. Recent numerical simulations (e.g., Haiman et al. 1996; Abel et al. 1998; 
Bromm et al. 1999) have shown that when the first collapsed objects with masses of 10 5-8 Mq are virialized, H2 relative 
abundance rises from its initial value of 10~ 6 to x^ 2 ~ 10~ 3 — 10 -4 and H 2 molecules then cool the gas to a temperature 
of T <~ 300 — 500 K. Consequently a thin baryonic disk forms, where the baryon density overwhelms that of the extended 
virialized dark halo (e.g. Umemura 1993). Then we have a relation such as 2nGpdH 2 ~ c^, where p&, H^, and c s denote 
the characteristic density, thickness of the disk, and sound speed, respectively. The density of the cooled disk is then 
estimated as 



Pd _ G Ml 



{j, 2nkT i? d 4 

» « * o&r m* & l d^r m"'^° • 



(A3) 

2 / 3 / O. \ 4 /3 



where p is the mean molecular weight, i?<j is the disk radius, and a(= Rd/r v i T ) is the disk radius in units of r v ; r . 
Furthermore, an overdense region acquires angular momentum through tidal spin-up by their surrounding fluctuations. 
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Then the overdense region cannot collapse into a disk smaller than the centrifugal barrier. The radius of the centrifugal 
barrier is given by 



r barr .0.01r vir ^j ^— j (A5) 

where A is a dimensionless spin parameter A = J|£ , | 1 / 2 G _1 M -5 / 2 (Sasaki & Umemura 1996), which is peaked around 
A = 0.05 (Heavens & Peacock 1988). The symbols J and E denote the total angular momentum and energy, respectively. 
The radius fbarr gives a lower bound for the size of the first collapsed pancake, which leads to a condition as a ;> 0.01. 

According to the linear theory of the disk in hydrostatic equilibrium, the disk is most unstable to the perturbation with 
the wavelength of A max = AttHa (e.g., Larson 1985). When the most unstable perturbation grows in the disk, the disk 
fragments into filamentary clouds rather than spherical clouds. Then, the separation of the filaments is equal to A max . 

If all the mass within one wavelength A max collapses into one filament of line mass I, then the mass of the filament is 
given as 

- ,3 x ^ ( £) (i±^)- (^) (^f ^y\- r . m 

The number of filaments formed from the parent disk is then estimated as 

*■ ~- 12 (£)"' (HH W (JkT (£f « ■ 
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Table 1 

REACTION RATE COEFFICIENTS 



Reactions 



Rate Cocfncicnts a 



Reference 



(HI) 
(H2) 
(H3) 
(H4) 

(H5) 

(H6) 

(H7) 
(H8) 

(H9) 

(H10) 
(Hll) 
(H12) 
(H13) 
(H14) 

(H15) 

(H16) 

(H17) 
(H18) 
(H19) 
(H20) 
(Hel) 
(Hc2) 



(Hc3) 



(He4) 

(He5) 
(He6) 



H + e- 
H+ + e 
H + H - 
H + e- 

H + H~ 

H + H+ 

H+ + H 
H 2 +H- 



H+ + 2e 
-+K + hv 
* H + H+ + e 
H~ + hv 



H 2 + e 



>H 2 
3H 



- hv 
H+ 



H 2 + H+ -» H+ + H 



H 2 

H 2 

H 2 

H 

H 

H 



e — 
e — 
H 2 
- e - 
H 



H + H- 

2H + e 
-» 2H + H 2 
>H + 2e 
-> 2H + e 



H+ -» 2H 



H- + H+ -» H+ + i 



H 2 + 
3H 
2H 
He- 



H 2 



-e-> 2H 
- H" -» H - 
^H 2 + H 

H 2 — » 2H 2 

e -> He+ + 2e 



He+ + e ^ He++ + 2e 



Hc+ + e -> He + hv 



He++ + e 

He + H+ 
He+ +H 



He 4 



hv 



He+ +H 
He + H+ 



fcff3 = 

km = 



Cff7 



5.85 x 10- n T - 5 exp(-157809.1/T)[l+ (T/10 5 ) 7 ]- 1 
8.40 x lO^T-^T/lO 3 )- - 2 !! + (T/IO 6 ) ' 7 ] -1 
1.7xl0- 4 fci 

1.4 x lCT 18 T - 928 cxp(-T/16200) 
1.5 x 10~ 9 T < 300K 

4.0 x 10- 9 T-° 17 T > 300K 
dex[-19.38- 1.523 log T 

+ 1.118(logT) 2 - 0.1269(logT) 3 ] 
6.4 x 10~ 10 



km [see eqs. (B6) — (B8) in reference] 



km = 



3.0 x 10~ 10 cxp (-21050/T) T < 10 4 K 
1.5 x 10~ 10 cxp (-14000/T) T > 10 4 K 
k H w = 2.7 x 10- 8 T- 127 cxp(-43000/T) 
k Hn = 4.38 x 10- 10 T - 35 exp(-102000/T) 
fc/fi2 [see eqs. (B6) — (B8) in reference] 



kmi 



k 



H15 = 

fcffl6 = 
fcffl7 = 

kms = 
knw = 
kma = 

fcffel = 
kRe2 = 



4.0 x 10- 12 cxp(-43000/T) 
5.3 x 10- 20 T 2 - 17 cxp(-8750/T) 
5.7 x 10- 6 T-°- 5 + 6.3 x 10~ 8 

-9.2 x 10- n T ' 5 + 4.4 x 10~ 13 j 
6.9 x i(r 9 T- - 35 T < 8000K 
9.6 x 10- 7 T-°- 9 T > 8000K 
2.0 x 10- 7 T-°- 5 
5.0 x 10- 6 T-°- 5 
5.5 x lO- 29 ^ 1 
kHw/8 

2.38 x 10- 11 r-°- 5 exp(-285335.4/T)[l- 
5.68 x 10- 12 T°- 5 cxp(-631515.0/T)[l 

3.294 x 10~ n 



f (T/10 5 ) - 5 ]" 1 
(T/10 5 ) 05 ]- 1 



kHe3 = 



\J 15 T 54 ( 1 + \J 15*54 ) 



0.309 



{} + \J 3.676x10^) 



1.691 



1.891 X 10 



-10 



KHe4 

kHe5 
kHe6 



T 

2.774X10 6 



^9.37 {} + ^9.37) 

) 



0.2476 



1.7524' 



+ 

4.0 x i(r 37 T 4 - 74 
3.7 x IQ-^T 2 - 06 [1 + 9.9 cxp (-3^0)] 



C92 
C92 
PSS83 
GP98 

GP98 

GP98 

GP98 
SK87 

GP98 

GP98 
SK87 
SK87 
SK87 
SK87 

GP98 

GP98 

GP98 

SK87 

PSS83 

PSS83 

C92 

C92 



VF96 



VF96 

GP98 
GP98 



a The units of rate coefficients are taken to be cm 3 s 1 for two-body reactions and cm 6 s 1 for three-body reactions. 



Note.— C92: Cen (1992); PSS83: Palla et al. (1983); GP98: Galli & Palla (1998); SK87: Shapiro & Kang (1987); 
VF96: Vcrncr & Ferland (1996) 
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Table 2 
MODEL PARAMETERS 



Model 


/ 


n c (cm J ) 


T (K) 


Xn 2 


Ala 


1.5 


10 


400 


10~ 4 


Alb 


1.5 


10 


300 


10~ 3 


A2a 


1.5 


10 2 


400 


10~ 4 


A2b 


1.5 


10 2 


300 


10~ 3 


A3a 


1.5 


10 3 


400 


10~ 4 


A3b 


1.5 


10 3 


300 


10~ 3 


A4a 


1.5 


10 4 


400 


10~ 4 


A4b 


1.5 


10 4 


300 


10~ 3 


A5a 


1.5 


10 5 


400 


10~ 4 


A5b 


1.5 


10 5 


300 


10~ 3 


A6a 


1.5 


10 6 


400 


10~ 4 


A6b 


1.5 


10 6 


300 


10~ 3 


Bla 


2 


10 


400 


10~ 4 


Bib 


2 


10 


300 


10~ 3 


B2a 


2 


10 2 


400 


10~ 4 


B2b 


2 


10 2 


300 


10~ 3 


B3a 


2 


10 3 


400 


10~ 4 


B3b 


2 


10 3 


300 


10~ 3 


B4a 


2 


10 4 


400 


io~ 4 


B4b 


2 


10 4 


300 


10~ 3 


B5a 


2 


10 5 


400 


10~ 4 


B5b 


2 


10 5 


300 


10~ 3 


B6a 


2 


10 6 


400 


io~ 4 


B6b 


2 


10 6 


300 


10~ 3 


Cla 


4 


10 


400 


io~ 4 


Clb 


4 


10 


300 


10~ 3 


C2a 


4 


10 2 


400 


io~ 4 


C2b 


4 


10 2 


300 


IO" 3 


C3a 


4 


10 3 


400 


io~ 4 


C3b 


4 


10 3 


300 


IO" 3 


C4a 


4 


10 4 


400 


io~ 4 


C4b 


4 


10 4 


300 


10~ 3 


C5a 


4 


10 5 


400 


io~ 4 


C5b 


4 


10 5 


300 


IO" 3 


C6a 


4 


10 6 


400 


io~ 4 


C6b 


4 


10 6 


300 


10~ 3 
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n (cm -3 ) n (cm 3 ) 

Fig. 1. — (a) Evolution of the temperature (T c ) and the H2 relative abundance (xh 2 ,c) at the center for model Ala with no density 
fluctuations. The model has the initial parameters of (no, To, /) = (10 cm -3 , 400 K, 1.5). For comparison, we also show the evolutional paths 
of the model in which the three-body reactions are neglected (dotted lines). The temperature stays nearly constant at 300 — 500 K over the 
seventh order of magnitude in density. When the central density reaches n c ~ 10 s cm -3 , the three-body reactions of H2 formation become 
dominant and almost all the hydrogen atoms become into H2 molecules. Thus, the H2 abundance steeply rises around n c ~ 10 8 cm -3 . When 
the density reaches n c ~ 10 12 cm -3 , the cloud becomes optically thick to the H2 lines, (b) Evolution of the contraction time (solid line), the 
cooling time (dashed line), and the fragmentation time (dashed- dotted line). For comparison, we also showed the evolution of the contraction 
time for the model in which the three-body reactions are neglected (dotted lines). The contraction becomes slower after the central density 
reaches the critical density of H2 (n C r ~ 10 3 cm -3 ). When the density reaches n c ~ 10 12 cm -3 , the cloud becomes optically thick to the H2 
lines, and the contraction essentially stops at that stage. For the model without the three-body reactions, the contraction decelerates at the 
stage at which n c ~ 10 9 cm -3 because of less effective H2 cooling. 




n (cm -3 ) n (cm 3 ) 



Fig. 2. — Same as Figure 1 but for model C6a. The model has the initial parameters of (no, To,/) = (10 6 cm - 3 , 400 K, 4). The overall 
evolution is qualitatively similar to that of model Ala. In this model, the contraction time does not exceed the fragmentation time until the 
cloud becomes optically thick to the H2 lines and the contraction essentially stops. 
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Fig. 3. — Density and velocity distributions in the r — z plane for model A4a at four different stages: (a) t = (n c = 1.0 X 10 4 cm 3 , 
initial state), (b) t = 5.9 X 10 5 yr (n c = 9.9 X 10 5 cm" 3 ), (c) t = 8.1 X 10 5 yr (n c = 5.1 X 10 7 cm" 3 ), and (d) t = 8.2 X 10 5 yr [n c = 7.5 X 10 7 
cm -3 ). At the initial state, scale-invariant density fluctuations with an amplitude of S = 0.1 are added. When the density reaches 10 7 cm" 3 , 
unstable density fluctuations begin to grow nonlinearly because the contraction time becomes longer than the fragmentation time. The mean 
mass and separation of the clumps are estimated as ~ 250Mq and 0.24 pc at the final stage of the computation, respectively. 
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Fig. 4. — Evolution of (a) the density and (b) velocity distributions in the z-axis for model A4a with scale-invariant density fluctuations 
of S = 0.1. They are depicted at the stages at which t = (solid lines), t = 5.9 X 10 s yr (long dashed lines), t = 8.1 X 10 5 yr (short dashed 
lines), and t = 8.2 X 10 5 yr (solid lines). When the density reaches 10 cm -3 , the density fluctuations begin to grow nonlinearly, and dense 
clumps form. In the clumps, the contraction is accelerated, and the temperature accordingly rises. 
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Fig. 5. — Evolution of (a) the density and (b) velocity distributions in the z-axis for model C6a with scale-invariant density fluctuations of 
5 = 0.1. They are depicted at the stages at which t = 5.458 X 10 4 yr (long dashed lines), t = 5.464 X 10 4 yr (dotted lines), and t = 5.47 X 10 4 
yr (solid lines). When the density reaches ~4x 10 12 cm -3 , the radial contraction essentially stops; then, the density fluctuations begin to 
grow nonlinearly, and dense clumps form. 
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Fig. 6. — Dependence of fragment mass on the initial central density of a filament no and a parameter / which corresponds to the initial 
line mass shown by the right ordinate. Scale-invariant (v = —1) density fluctuations with an amplitude of 8 = 0.1 are assumed. Solid lines 
denote the contour lines of the averaged fragment mass. Dashed lines show the line masses of the filaments formed by cosmological collapse 
and fragmentation (see text and the Appendix). The numbers attached to the dashed lines denote the line masses normalized with the 
equilibrium value (e.g., l/l cq = 1, 1.5, and 2). When the line mass of the filament is less than l/l C q = 1, the cloud expands because the 
pressure force is greater than the gravitational force (light-gray region). The averaged fragment mass depends on the initial parameters. For 
larger / and/or higher initial density, the fragment mass is lower. The maximum and minimum masses are estimated as 10 Mq and 1 ~ 2 
Mq, respectively. These two masses are related to the microphysics of H2 molecules. The former corresponds to the Jeans mass at the stage 
at which the density reaches the critical density of H2 molecules, while the latter corresponds to the Jeans mass at the stage at which the 
cloud becomes opaque to the H2 lines. There is a steep boundary in the fragment mass at no ~ 10 5-6 cm -3 and / > 3. For the models with 

no > 10 5 cm -3 , the fragment masses take their minimum at 1 ~ 2 Mq, while for the models with no < 10 s cm -3 , they are greater than 

~ 10 2 Mq. This sensitivity in the fragment mass comes from the rapid increase in the H2 abundance due to the three-body reactions. 
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Fig. 7. — Same as Figure 6 but for 8 = 0.01. This figure shows that the dependence on the fluctuation amplitude is quite weak and the 
fragment mass is not basically altered. 
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Fig. 8. — Schematic initial mass function of Population III stars. Solid and dashed lines are the IMFs for 3cr and ltr density fluctuations, 
respectively. The IMF of Population III stars is likely to be bimodal and approximated by a superposition of two power-law like components 
with two different peaks of m p i « 1 — 2Mq and m p 2 ~ a few X 10 — 1O 2 M0. The relative height of the first peak descends with time compared 
to the second peak. According to a recent theory of stellar evolution, stars with masses between 1 — 2Mq and 8Mq are likely to evolve 
to white dwarfs that may reside in galactic halos as baryonic dark matter. Stars with masses between 8Mq and 35Mq probably evolve to 
supernovae and eject heavy elements in the intergalactic medium. Stars with masses between 35Mq and 10 2 Mq and greater than 25OM0 are 
likely to collapse into black holes and may be responsible for baryonic dark matter. Stars with masses between 1O 2 M0 and 25OM0 probably 
explode as supernovae and inject heavy elements into the intergalactic medium. 



